Matter-wave dark solitons in a double-well potential 
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We study stability of the first excited state of quasi-one-dimensional Bose-Einstein condensates in 
a double- well potential, which is called "7r-state" . The density notch in the 7r-state can be regarded 
as a standing dark soliton. From the excitation spectrum, we determine the critical barrier height, 
above which the 7r-state is dynamically unstable. We find that the critical barrier height decreases 
monotonically as the number of condensate atoms increases. We also simulate the dynamics of the n- 
state by solving the time-dependent Gross-Pitaevskii equation. We show that due to the dynamical 
instability the dark soliton starts to move away from the trap center and exhibits a large-amplitude 
oscillation. 

PACS numbers: 03.75.Hh, 03.75.Kk, 03.75.Lm, 05.30.Jp 



I. INTRODUCTION 



Since the first discovery of Bose-Einstein condensates 
(BECs) of alkali atoms in 1995 [HQ, macroscopic quan- 
tum phenomena, which stem from the coherence of the 
many-body wavefunctions describing BECs, have been 
studied extensively. Systems of a BEC in a double-well 
potential arc well suited for exploring various phenomena 
associated with the coherent nature, and vigorous studies 
in such systems in recent years have led to the observa- 
tion of intriguing macroscopic quantum phenomena, such 
as the matter- wave interference [1| and the Josephson ef- 
fects H] . 

The coherent nature of Bose-condensed atomic gases 
combined with their diluteness allows for the formulation 
of the BEC systems based on a non-linear Schrodingcr 
equation, the so-called Gross-Pitaevskii (GP) equa- 
tion The static solutions of the GP equation de- 
scribe a BEC not only in the ground state but also in 
the macroscopically excited states, such as the station- 
ary quantized vortices @ and the standing dark soli- 
tons. In particular, the dark solitons in BECs, which 
share fundamental properties with solitons in the nonlin- 
ear optics and the fluid mechanics, have attracted much 
attention 0, d, @ • Most recently, Becker et al. observed 
experimentally dark solitons with long lifetime up to sev- 
eral seconds and opened up newpossibility of detailed 
studies of the soliton dynamics j9[. Previous theoreti- 
cal work [l(J [ll[ on stability of the dark soliton states 
in harmonic potentials with axial symmetry found that 
the dark soliton is stable when the radial confinement 
is sufficiently strong, i.e. when the BEC is quasi-one- 
dimensional (ID). On the other hand, in 3D systems a 
dynamical instability, which is called the snake instabil- 
ity, causes the breakdown of the dark soliton, resulting 
in the formation of vortex rings 0, HH • 

The first excited state of a BEC in a double-well po- 
tential has a density notch and a phase kink at the center 
of the trap (see the dashed line in Fig. [1]) . This state is 
called 7r-state because the phase difference between the 
left and right condensates is ir. In the sufficiently large 
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FIG. 1: (Color online) The condensate wavefunctions for the 
ground state (solid line) and the 7r-state (dashed line) of 
BECs in a double-well potential. We set the potential barrier 
Vo/(?UjJz) = 10.0. We consider condensates of 23 Na atoms, 
the mass is 23 amu (atomic mass units) and the atom num- 
ber N = 2000. The values of parameters are as follows: the 
harmonic-oscillator frequency is u) z /(2n) = 10 Hz, the fre- 
quency of the radial confinement is ui ±/(2n) = 250 Hz and 
the s-wave scattering length is a s — 3 nm. 



condensate, the size of the density notch becomes of the 
order of the healing length. In this regime, the density 
notch can be regarded as a dark soliton. Stability of the 
7r-state of quasi-lD BECs has been previously studied 
only in two extreme limits; it is dynamically stable when 
the potential barrier is absent [lj| , while it is dynamically 
unstable when the potential barrier height is much larger 
than the chemical potential, i.e. when the two-mode ap- 
proximation is valid [l2j|. Although the transition from 
the dynamically stable 7r-state to the dynamically un- 
stable 7r-state is expected to occur, no work to date has 
studied the intermediate transition region. In order to 
identify the transition, one needs linear stability analy- 
ses based on the Bogoliubov equations without using the 
two-mode approximation, which is not applicable in the 
transition region. 

In this paper, we study stability of the 7r-state at zero 
temperature in a wide range of the barrier height on the 
basis of the GP mean- field theory. First, we calculate the 
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excitation energies of the 7r-state by solving the Bogoli- 
ubov equations numerically. The critical barrier height, 
above which the 7r-state is dynamically unstable, is de- 
termined from the condition of complex excitation en- 
ergy. We will show that the critical barrier height de- 
creases monotonically as the number of atoms increases. 
For comparison, we also calculate the lowest excitation- 
energy using the two-mode approximation. This calcula- 
tion will confirm that the two-mode approximation com- 
pletely fails to provide the critical barrier height except 
in the limit of small number of atoms. Secondly, we nu- 
merically simulate the dynamics of the 7r-state by solving 
the time-dependent GP equation. We will find that in 
the dynamically unstable region the condensate density 
notch corresponding to the dark soliton oscillates with a 
large amplitude around the center of the trap, while in 
the dynamically stable region the dark soliton remains 
stationary at the trap center during the time evolution. 

The present paper is organized as follows. In Sec. [Til 
we start with the model and formulation based on the 
GP mean-field theory. In Sec. IIII| we solve the time- 
independent GP equation and the Bogoliubov equations 
numerically, and we obtain the stability phase diagram. 
In Sec. |TVj we numerically solve the time-dependent GP 
equation to elucidate the decay process of the dynami- 
cally unstable condensates. In Sec. El we summarize our 
results. 



II. MODEL AND FORMULATION 

We consider a BEC at zero temperature confined in 
a double-well potential. The double-well potential is as- 
sumed to be composed of a harmonic trap with cylindri- 
cal symmetry and a Gaussian-shaped potential barrier, 
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where u>± (lo z ) is the frequency of the radial (axial) har- 
monic confinement, Vq is the barrier height and a is the 
width of the barrier. In typical experiments, a potential 
barrier is created by a blue-detuned laser beam and its 
height is controlled by the intensity of the laser beam. 
We assume that huj± is much larger than the chemical 
potential so that we can justify the ID treatment of this 
problem. Accordingly, equilibrium and noncquilibrium 
properties of BECs in a double-well potential are de- 
scribed by the ID time-dependent GP equation for the 
condensate wavefunction ^{z, t): 
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The coupling constant g in ID is given by g — 
2h 2 a s /ma 2 l , where a s is the s-wave scattering length and 
a± = (h/muj±) 1 ^ 2 is the harmonic oscillator length of the 
radial confinement. 

The static solution tf CM) = ^n{z)&^ t/h of Eq. © 
satisfies the time-independent GP equation, 
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where fi is the chemical potential. &o(z) satisfies the 
normalization condition, 



\$ (z)fdz = N, 



(5) 



where N is the number of condensate atoms. Considering 
small oscillations around the static solution <!>o(z), we 
write the condensate wavefunction as 
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Substituting Eq. (JSJ) into Eq. @ and linearizing it with 
respect to the quasiparticle amplitudes u{z) and v(z), we 
obtain the Bogoliubov equations 
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where we introduced the operator 
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The solutions of the Bogoliubov equations ([7]) determine 
the excitation energy e. 

We solve the time-independent GP equation (fj| and 
the Bogoliubov equations ([7]) using the method intro- 
duced in Ref. [l3j]. We employ the harmonic oscillator 
cigenfunctions <p\{z) as a basis set, where 
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They satisfy the orthonormality relation 

dzip* x (z)(fi X '(z) = S\\>- 
We expand the static solution as 

$oO) = \/iV^a A ^A(;z), 
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and insert Eq. (TT2]) into Eq. (TJ}. Multiplying Eq. Q by 
fx( z ) an< l making use of Eq. (|11[) . we obtain 
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FIG. 2: Stability phase diagram calculated from the Bogoli- 
ubov equations (0 (solid line) and the two-mode approxima- 
tion (dashed line). 



One can solve these nonlinear simultaneous equations for 
cl\ by the Newton-Raphson method. Since the conden- 
sate wavefunction of the 7r-state is an odd function, we 
take only the harmonic oscillator cigcnfunctions with odd 
parity in order to obtain the 7r-state. For instance, in 
Fig. [1] we show the ground state and the 7r-state for N = 
2000 and Vo/(tiu) z ) = 10.0. Subsequently, by substituting 
the condensate wavefunction $o( z ) m to the Bogoliubov 
equations ([7]) and expanding u and v in terms of the har- 
monic oscillator eigenfuntions as u{z) = u\(p\(z) and 
v ( z ) — J2\ v ^ l P^( z )j we obtain the simultaneous linear 
equations for u\ and v\. One can obtain the excitation 
energy e by diagonalizing the coefficient matrix of the 
linear equations. 

In the following sections, we study stability of the 7r- 
state in a double-well potential. As shown in Fig. [JJ the 
major difference between the ground state and the w- 
state is that the 7r-state has a notch with zero density 
and the phase jump at the center of the trap. The notch 
of the condensate density with the phase jump of it is 
regarded as a dark soliton, when the number of atoms is 
large (Thomas-Fermi limit [H ) . In the case of the ground- 
state BEC, the relation between the excitation energies 
of the ground state and the barrier heig ht in a double- 
well potential has been studied [3, Il5j . It was found 
that the lowest excitation-energy exhibits the crossover 
from the dipole mode to the Josephson plasma mode as 
the barrier height increases 



III. STABILITY ANALYSIS 

In this section, we solve the time-independent GP 
equation ([¥]) and the Bogoliubov equations ([7]) numer- 
ically and obtain the stability phase diagram. The sta- 
bility of BECs can be studied by calculating the exci- 
tation energy. The appearance of excitations with com- 
plex energies signals the dynamical instability whereby 
the amplitude of the condensate fluctuation grows expo- 
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FIG. 3: The lowest or imaginary eigenvalue of the Bogoliubov 
equations ((7} as a function of Vo/{hui z ) with the fixed atom 
number N = 2000. 
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FIG. 4: The closed circles represent the excitation energy for 
N — 2000 as a function of Vo/(Hu) z ). etma is also plotted by 
the opened circles. 



nentially in time. On the other hand, the appearance of 
excitations with negative energies signals the Landau (en- 
ergetic) instability, which means that the solution $o of 
the time-independent GP equation does not correspond 
to a local minimum in the energy landscape. Neverthe- 
less, the Landau instability cannot destabilize the system 
in the limit of zero temperature because of the absence of 
dissipative processes [H], EH E3 ■ In other words, the sys- 
tem is dynamically stable under the condition Re[e] < 
and Im[e] = 0. 

Fig-IHshows the stability phase diagram in the (N, Vq)- 
plane. The solid line represents the critical barrier height 
above which Im[e] ^ so that the system is dynamically 
unstable. One sees that the critical barrier height mono- 
tonically decreases as N increases. 

From the numerical solution of Eq. ([JJ), we find that 
there is only one eigenstate with a negative excitation 
energy in the dynamically stable re gion (as also found 
in a special case of Vq = in Ref. |10[)- This excita- 
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FIG. 5: The closed circles represent the excitation energy for 
TV = 10 as a function of Vo/(hLj z ). etma is also plotted by the 
opened circles. 



tion corresponds to the oscillation of the density notch 
around the trap center. The reason for the negative en- 
ergy for this excitation can be easily understood for the 
case of Vb = 0. In this case, it is clear that the potential 
energy for the density notch is the highest at the cen- 
ter of the trap. Thus, the deviation of the density notch 
from the potential maximum lowers the energy, resulting 
in the negative excitation energy. As Vq increases for a 
fixed value of N, the absolute value of this negative exci- 
tation energy decreases because the potential energy for 
the density notch also decreases. As Vq increases from 
the critical point, the imaginary part of the excitation 
energy develops and takes the maximum value at a cer- 
tain point. With increasing barrier height Vq further, the 
imaginary part of the excitation energy decreases; it ap- 
proaches zero asymptotically in the limit Vq — > oo. In 
this limit, the two condensates separated by the poten- 
tial barrier are completely independent of each other and 
thus the system becomes dynamically stable again. 

Here we discuss the case of N = 2000, where the mean- 
field interaction is so strong that the density notch of the 
condensate at the center of the trap can be regarded as a 
dark soliton. The excitation energy of the mode relevant 
to the stability of the 7r-state is shown in Figs. [3] and 
31 The trends mentioned above are clearly seen in these 
figures. The critical barrier height is Vo/(huj z ) ~ 0.24, 
while the barrier height that gives the maximum value of 
Im[e] is V /(hu} z ) -6.0. 

We next calculate the excitation energy within the two- 
mode approximation. This approximation is well known 
to provide analytical insights into BECs in a double-well 
potential when the barrier height is sufficiently large com- 
pared to the chemical potential [H, [3, EH ■ The analyt- 
ical expression of the excitation energy can be written 
as 



which is easily derived from the Joscphson Hamilto- 
nian [13 ]. The coupling energy K and the capacitive 
energy U arc defined as 
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where E. 



sm(an 



is the energy of the condensate of the sym- 
metric (antisymmetric) state, and fi sra is the chemical 
potential of the symmetric state. The energy of the con- 
densate is given by 
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The dashed line in Fig. [5] represents the transition 
point between the dynamically stable and unstable re- 
gions determined from the condition Im[e tma ] =/= 0. In 
Fig. O £tma gives almost the same value as the excitation 
energy e determined from the Bogoliubov equations ([7J 
even near the critical barrier height, since the transition 
point is located in the region where the barrier height 
is larger than the chemical potential. In contrast, it is 
clearly seen in Fig. 5] that e t ma for A^=2000 completely 
disagrees with e in the vicinity of the transition point, 
where Vq -C /i- Thus, the two-mode approximation fails 
to provide a correct phase boundary in the large N re- 
gion. 



IV. DYNAMICS OF THE CONDENSATE WAVE 
FUNCTION 

In this section, we solve the time-dependent GP equa- 
tion (J2J) numerically to elucidate the decay process of the 
dynamically unstable condensates. We use the split op- 
erator fast-Fourier transform method. As in Fig. [3l we 
set the number of atoms as N = 2000. In order to trigger 
the instability, we need to add a small random noise to 
the initial state as 



(17) 



The amplitude of the random noise is chosen to be 
J dz|<5'I'| 2 — N x 10~ 4 . In actual experiments, there al- 
ways exists such a small noise due to thermal and quan- 
tum fluctuations of the condensates or the vibration of 
the external field potential. 

In Figs. [6] and [7j we plot the time evolution of the den- 
sity ^(z,^)! 2 and phase <fr(z,t) of the condensate. The 
phase (j>(z,t) is an important quantity in understanding 
the soliton dynamics since the phase jump at the density 
notch A(j) is related to the velocity of the gray soliton v 
by 
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e t „ la - y/2K(2K-UN), 
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where c is the sound velocity Moreover, we show the 
position of the dark soliton as a function of time in Fig. [5] 
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FIG. 6: (Color online) Time evolution of (a) the density and 
(b) the phase of the condensate wavefunction obtained by 
solving the time-dependent GP equation starting with the 7r- 
state. The barrier height is Vo/(fiui z )—0, where the system 
is dynamically stable. The dimensionless time is defined by 

t = tliJz- 
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FIG. 7: (Color online) Time evolution of (a) the density and 
(b) the phase of the condensate wavefunction in the 7r-state 
obtained by solving the time-dependent GP equation in the 
dynamically unstable region Vb/(?iWz)=6.0. The dimension- 
less time is defined by t = tu) z . 
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FIG. 8: (Color online) The kink position z of the gray soliton 
as a function of the dimensionless time t = tu z . The cases 
of Vo/(hui z )=0, 0.1 and 0.2 (in the dynamically stable region) 
correspond to the solid line. The dashed, dotted and dashed- 
dotted lines represent the kink positions for Vo/(hu> z )—2.0, 
Vo/(hio z )—4.0 and Vo/(huj z )—6.0 (in the dynamically unstable 
region), respectively. 



As seen in Fig. [5] and the solid line in Fig. [51 when 
Vo/(huj z ) = 0, the density notch remains stationary at 
the trap center during the time evolution and the phase 
jump remains 7T. Thus, the system is dynamically stable. 

When the barrier height is above the critical value as 
shown in Fig. [7] (Vo/(Hlj z ) = 6.0) and the dashed-dotted 
line in Fig.[SJ the density notch starts to move away from 
the center of the trap after a certain time t ~ (Im[e/?i]) _1 
owing to the dynamical instability. The density notch 
then develops a large-amplitude oscillation. One can see 
that the notch moves with larger speed as the absolute 
value of the phase jump becomes smaller, consistent with 
Eq. (QID. 

Since the condensate dynamics in the dynamically un- 
stable region is distinctively different from that in the sta- 
ble region as discussed above, one expects that the transi- 
tion to the dynamically unstable region can be identified 
experimentally. We suggest the following procedure to 
observe the transition in experiments: at first one could 
engineer the 7r-state in a double- well trap by means of the 
phase imprinting techniques 0, [l9| . Subsequently, one 
could tune the barrier height by controlling the intensity 
of the blue-detuned laser beam. When the barrier height 
is so small that the system is in the stable region, one 
should observe a dark soliton localized at a trap center. 
Setting the barrier height to be large enough to exceed 
the critical value, one should observe a large-amplitude 
oscillation of a dark soliton. Thus, the transition should 
be identified by measuring the time evolution of the con- 
densate density profile. 



V. CONCLUSION 

In conclusion, we have studied stability of the first ex- 
cited state of Bose-Einstcin condensates in a double-well 
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potential on the basis of the Gross-Pitaevskii mean-field 
theory. From the excitation energies obtained from the 
Bogoliubov equations, we determined the critical barrier 
height above which the system is dynamically unstable. 
The stability phase diagram was obtained against the 
barrier height Vq and the number of condensate atoms 
N. We found that the critical barrier height monotoni- 
cally decreases as N increases. 

Solving the time-dependent GP equation numerically, 
we also studied the time evolution of the condensate in 
the 7r-state. Our simulation results show that the con- 
densate density notch in the dynamically unstable region 
develops a large-amplitude oscillation; this behavior is 



significantly different from that in the dynamically sta- 
ble region. 
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